!**********************************
! -----------------------------------------------------------------------
! File:             elliptic_solvers.f90
! Purpose:          routines related to the solution of elliptic solvers
! Contains:
! Revision History: Ver. 1.0 May. 2007 Alberto Scotti
! Added tridiagonal solver May 2008 Alberto Scotti
! -----------------------------------------------------------------------





 

!********************
MODULE Elliptic_Solvers
!********************
! USER definitions of data structures associated with each node
USE ElliDef
USE metric
USE MAtrix_Oper

IMPLICIT NONE
SAVE
PRIVATE

INTEGER :: maxfct,mnum,mtype,phase,msglvl,nrhs,idum,error

REAL (KIND=Precision) :: ddum

INTEGER, PARAMETER :: BEFORE_POISSON=10,AFTER_POISSON=11,ErrOk=0, EXTR_BC=0

INTEGER, PARAMETER, PUBLIC :: EXTR_BC2d=10,EXTR_BC3d=0
INTEGER, PUBLIC, PARAMETER :: CHECK_BC=1000, &
     CHECK_COEFFICIENTS = 1001,SOLVE=10000,INIT=10001,FINISH=10002 !these are defined in
                                          ! nodeInfoGlobal and recycled here
INTEGER, PUBLIC, PARAMETER :: IT=-1,IT_DBCG=-2,IT_BCG=-3,IT_CG=-5,&
     IT_CGNR=-7,IT_FOM=-9,IT_TFQMR=-10,IT_BCGSTAB=-11,IT_DQGMRES=-12,IT_GMRES=-13,IT_FGMRES=-14,IT_SD=-15,PARDISO_SOLVER=100,     CG_SOLVER=101,DCG_SOLVER=102

PUBLIC project,take_div,check,dcg_solve,POISSON_SOLVERS,CREATE_OPERATORS,FREEUP_OPERATORS,POISSON_SOLVERS_THIN,CREATE_OPERATORS_THIN,APPLY_BC,hodge_thin,hodge_thick,set_solver_parms,get_solver_parms

CONTAINS 

SUBROUTINE SET_SOLVER_PARMS(Info,which,what)
TYPE (Elli) :: Info
INTEGER, OPTIONAL, INTENT (IN) :: which,what
! used to modify parameters of the solver from defaults
IF(.NOT.PRESENT(which)) THEN
   PRINT*,'set_solver_parms invoked without specifying which parameter to modify'
   STOP
ENDIF
IF(.NOT.PRESENT(what)) THEN
   PRINT*,'set_solver_parms invoked without specifying what value to use'
   STOP
END IF
IF(which<0.OR.which>64) THEN
   PRINT*,'which in set_solver_parms exceeds bounds'
   STOP
END IF
Info%it%ipar_pd(which)=what
END SUBROUTINE SET_SOLVER_PARMS

SUBROUTINE GET_SOLVER_PARMS(Info,which,what)
TYPE (Elli) :: Info
INTEGER, OPTIONAl, INTENT (IN) :: which
INTEGER, OPTIONAL, INTENT (OUT) :: what
! used to query the value of a parameter of the solver
IF(.NOT.PRESENT(which)) THEN
   PRINT*,'get_solver_parms invoked without specifying which parameter to modify'
   STOP
ENDIF
IF(.NOT.PRESENT(what)) THEN
   PRINT*,'get_solver_parms invoked without specifying what value to use'
   STOP
END IF
IF(which<0.OR.which>64) THEN
   PRINT*,'which in get_solver_parms exceeds bounds'
   STOP
END IF
what=Info%it%ipar_pd(which)

END SUBROUTINE GET_SOLVER_PARMS

SUBROUTINE setparm(Info,solver) ! eventually move this to problem...
TYPE (Elli) :: Info
INTEGER, EXTERNAL :: omp_get_max_threads
INTEGER, INTENT (IN) :: solver
TYPE (CRS), POINTER ::  A=>NULL()
SELECT CASE (solver)

CASE (PARDISO_SOLVER)
ALLOCATE(Info%it%ipar_pd(64))
! 0.  
! allocation of resources
   Info%it%ipar_pd=0

   Info%it%ipar_pd(1)=1 ! no solver default
   Info%it%ipar_pd(2)= 2 ! fill-in reordering from METIS
   Info%it%ipar_pd(3) = omp_get_max_threads() !numbers of processors, value of OMP_NUM_THREADS
   Info%it%ipar_pd(4) = 0 ! no iterative-direct algorithm
   Info%it%ipar_pd(5) = 0 ! no/yes user fill-in reducing permutation
   Info%it%ipar_pd(6) = 0 ! =0 solution on the first n compoments of x
   Info%it%ipar_pd(7) = 16 ! default logical fortran unit number for output
   Info%it%ipar_pd(8) = 2 ! numbers of iterative refinement steps
   Info%it%ipar_pd(9) = 0 ! not in use
   Info%it%ipar_pd(10) = 10 ! perturbe the pivot elements with 1E-13
   Info%it%ipar_pd(11) = 0 ! use nonsymmetric permutation and scaling MPS
   Info%it%ipar_pd(12) = 0 ! not in use
   Info%it%ipar_pd(13) = 0 ! not in use
   Info%it%ipar_pd(14) = 0 ! Output: number of perturbed pivots
   Info%it%ipar_pd(15) = 0 ! not in use
   Info%it%ipar_pd(16) = 0 ! not in use
   Info%it%ipar_pd(17) = 0 ! not in use
   Info%it%ipar_pd(18) = -1 ! Output: number of nonzeros in the factor LU
   Info%it%ipar_pd(19) = -1 ! Output: Mflops for LU factorization
   Info%it%ipar_pd(20) = 0 ! Output: Numbers of CG Iterations
   error = 0 ! initialize error flag
   msglvl = 0 ! don't print statistical information
   mtype = -2 ! symmetric
   
   
   nrhs=1;maxfct=1;mnum=1
   CASE (IT)
      ALLOCATE(Info%it%ipar_it(16),Info%it%fpar(16))
      Info%it%ipar_it=0
      Info%it%fpar=zero

      Info%it%ipar_it(2)=1 ! use left preconditioner
      Info%it%ipar_it(3)=1
      Info%it%ipar_it(4)=-999 ! to be filled later
      Info%it%ipar_it(5)=15 ! size of Krylov subspace
      Info%it%ipar_it(6)=15 ! maxiteration
      Info%it%fpar(1)=1.0d-3
      Info%it%fpar(2)=1.d-12

      Info%it%droptol=1.d-3
      Info%it%lfil=10
CASE DEFAULT

   PRINT*,'setparm called with unknown request'
   STOP
END SELECT

END SUBROUTINE setparm

SUBROUTINE pardiso_init(Info,oper)
! Allocates memory and computes the LU decomposition of the matrix.
! to be called every time a new grid is generated

TYPE (Elli) :: Info
INTEGER (KIND=INT_PTR_KIND()), POINTER, DIMENSION (:) :: pt
INTEGER, POINTER, DIMENSION (:) :: columns,rowIndex
REAL (KIND=Precision), POINTER, DIMENSION(:) :: rValues
INTEGER :: i,j,k,l,m
INTEGER, OPTIONAL :: oper
INTEGER :: N,OP
! set the parameters of the solver

N=Product(Info%mX(1:Info%nDim))
! allocate memory pointers
OP=CSR_LAPL
!!$msglvl=1
IF(PRESENT(oper)) THEN
   OP=oper
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   Info%it%ipar_pd(5)=1 ! if 1, uncomment lines below to
                        ! define perm modify value from standard
!!$   msglvl=0 ! one for diagnostic from LU factor
   ! with this permutation, the precon matrix (assumed it has been calculated
   ! setting to zero all but the z-derivatives) becomes a tridiagonal matrix
   ! This eliminates the expensive step of having to calculate the 
   ! permutation. For a general preconditioner DO NOT do this...
   ALLOCATE(Info%it%perm(PRODUCT(Info%mX(1:3))))
   DO k=1,Info%mX(3)
      DO j=1,Info%mX(2)
         DO i=1,Info%mX(1)
            l=(k-1)*Info%mX(1)*Info%mX(2)+(j-1)*Info%mX(1)+i
            Info%it%perm(l)=(i-1)*Info%mX(3)*Info%mX(2)+(j-1)*Info%mX(3)+k
         ENDDO
      ENDDO
   ENDDO
ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF(PRESENT(oper)) THEN ! user ptprcond
   IF(ASSOCIATED(Info%ptprcond)) THEN
      DEALLOCATE(Info%ptprcond);nullify(Info%ptprcond)
   ENDIF
   ALLOCATE(Info%ptprcond(64));
   pt=>Info%ptprcond
   pt=0

ELSE
   IF(ASSOCIATED(Info%pt)) THEN
      DEALLOCATE(Info%pt);nullify(Info%pt)
   ENDIF
   
   ALLOCATE(Info%pt(64));
   pt=>Info%pt
   pt=0
   
END IF
 IF(.NOT.ASSOCIATED(Info%G(OP)%crsptr%rValues)) THEN
    PRINT*,'sparse matrix not present'
    PRINT*,'pardiso_init exits now'
    STOP
 ENDIF

rValues=>Info%G(OP)%crsptr%rvalues;columns=>Info%G(OP)%crsptr%columns;
rowIndex=>Info%G(OP)%crsptr%rowIndex
! initialize pardiso
phase=12 ! analysis + factorization

   
   CALL PARDISO(pt,maxfct,mnum,mtype,phase,N,rValues,rowIndex,columns, &
        Info%it%perm,nrhs,Info%it%ipar_pd,msglvl,ddum,ddum,error)

IF(error/=0) THEN
   PRINT*,'error',error,' occurred during pardiso init'
   PRINT*,'exiting now'
   STOP
ENDIF
msglvl=0;
nullify(pt,rvalues,columns,rowIndex)
END SUBROUTINE pardiso_init

SUBROUTINE pardiso_mopup(Info,oper)
! release all internal memory used for pardiso solver
! it deallocates Info%pt, but not the storage for the sparse
! matrix
TYPE (Elli) :: Info
INTEGER (KIND=INT_PTR_KIND()), POINTER, DIMENSION (:) :: pt
INTEGER, POINTER, DIMENSION (:) :: columns,rowIndex
REAL (KIND=Precision), POINTER, DIMENSION(:) :: rValues
INTEGER :: N,OP
INTEGER, OPTIONAL :: oper
N=Product(Info%mX(1:Info%nDim))
! allocate memory pointers

IF((.NOT.ASSOCIATED(Info%pt)).AND.(.NOT.ASSOCIATED(Info%ptprcond))) THEN
   PRINT*,'cannot execute pardiso mopup'
   PRINT*,'pointers not allocated',ASSOCIATED(Info%pt),ASSOCIATED(Info%ptprcond)
   PRINT*,'exiting now'
   STOP
ENDIF
OP=CSR_LAPL
IF(PRESENT(oper)) THEN
   OP=oper
   pt=>Info%ptprcond
ELSE
   pt=>Info%pt
ENDIF
 IF(.NOT.ASSOCIATED(Info%G(OP)%crsptr%rValues)) THEN
    PRINT*,'sparse matrix not present'
    PRINT*,'pardiso_init exits now'
    STOP
 ENDIF

rValues=>Info%G(OP)%crsptr%rvalues;columns=>Info%G(OP)%crsptr%columns;
rowIndex=>Info%G(OP)%crsptr%rowIndex
! mopup pardiso
phase=-1 ! release all internal memory
CALL PARDISO(pt,maxfct,mnum,mtype,phase,N,rValues,rowIndex,columns, &
     idum,nrhs,Info%it%ipar_pd,msglvl,ddum,ddum,error)
IF(ASSOCIATED(Info%it%perm)) DEALLOCATE(Info%it%perm)
IF(error/=0) THEN
   PRINT*,'error',error,' occurred during pardiso mopup'
   PRINT*,'exiting now'
   STOP
ENDIF

IF(PRESENT(oper)) THEN
   DEALLOCATE(Info%ptprcond);nullify(Info%ptprcond);
ELSE
   DEALLOCATE(Info%pt);nullify(Info%pt);
ENDIF
nullify(pt,rvalues,columns,rowIndex)

END SUBROUTINE pardiso_mopup

SUBROUTINE pardiso_solution(Info,oper,lhs,rhs)

! solves the ellyptic problem

TYPE (Elli) :: Info
INTEGER, OPTIONAL :: oper
REAL (KIND=Precision), DIMENSION (:), OPTIONAL :: lhs,rhs
INTEGER (KIND=INT_PTR_KIND()), POINTER, DIMENSION (:) :: pt,pt_temp
INTEGER, POINTER, DIMENSION (:) :: columns,rowIndex
REAL (KIND=Precision), POINTER, DIMENSION(:) :: rValues
REAL (KIND=Precision), POINTER, DIMENSION(:,:,:) :: p,div
INTEGER :: N,OP

OP=CSR_LAPL
IF(PRESENT(oper)) THEN
OP=oper
pt=>Info%ptprcond
ELSE
pt=>Info%pt
ENDIF
N=Product(Info%mX(1:Info%nDim))
IF(.NOT.ASSOCIATED(Info%pt).AND..NOT.ASSOCIATED(Info%ptprcond)) THEN
   PRINT*,'cannot execute pardiso solve'
   PRINT*,'pointers not allocated'
   PRINT*,'exiting now'
   STOP
ENDIF
   

 IF(.NOT.ASSOCIATED(Info%G(OP)%crsptr%rValues)) THEN
    PRINT*,'sparse matrix not present'
    PRINT*,'pardiso_solve exits now'
    STOP
 ENDIF

rValues=>Info%G(OP)%crsptr%rvalues;columns=>Info%G(OP)%crsptr%columns;
rowIndex=>Info%G(OP)%crsptr%rowIndex

phase=33 ! solution and iterative refinement
IF(PRESENT(lhs)) THEN
CALL PARDISO(pt,maxfct,mnum,mtype,phase,N,rValues,rowIndex,columns, &
     idum,nrhs,Info%it%ipar_pd,msglvl,rhs,lhs,error)

ELSE

div=>Info%div
p=>Info%p(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)) 

! solve

CALL PARDISO(pt,maxfct,mnum,mtype,phase,N,rValues,rowIndex,columns, &
     idum,nrhs,Info%it%ipar_pd,msglvl,div,p,error)
!!$      ! calculate and print residual
!!$      p=>adv_sym(Info%G(OP)%crsptr,p,-(/0,0,0/))
!!$      
!!$      PRINT 1000,sqrt(SUM((p-div)**2)/PRODUCT(Info%mX(1:3)))/&
!!$           SQRT(SUM(div**2)/PRODUCT(Info%mX(1:3)))
!!$      PRINT*,Info%it%ipar_pd(4)
!!$      DEALLOCATE(p)


ENDIF
IF(error/=0) THEN
   PRINT*,'error',error,' occurred during pardiso solution'
   IF(error==-4) THEN
      PRINT*,'value of iparm(20) is',Info%it%ipar_pd(20)
      PRINT*,'I will recalculate'
      p=>Info%p(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)) 
      p=zero
      ALLOCATE(pt_temp(64))
      pt_temp=0
      Info%it%ipar_pd(4)=0

      CALL PARDISO(pt_temp,maxfct,mnum,mtype,12,N,rValues,rowIndex,columns, &
     Info%it%perm,nrhs,Info%it%ipar_pd,0,div,p,error)
      CALL PARDISO(pt_temp,maxfct,mnum,mtype,33,N,rValues,rowIndex,columns, &
     Info%it%perm,nrhs,Info%it%ipar_pd,0,div,p,error)

      CALL PARDISO(pt_temp,maxfct,mnum,mtype,-1,N,rValues,rowIndex,columns, &
     Info%it%perm,nrhs,Info%it%ipar_pd,msglvl,div,p,error)
      DEALLOCATE(pt_temp)
            p=>adv_sym(Info%G(OP)%crsptr,p,-(/0,0,0/))
      
      PRINT 1001,sqrt(SUM((p-div)**2)/PRODUCT(Info%mX(1:3)))/&
           SQRT(SUM(div**2)/PRODUCT(Info%mX(1:3)))
      DEALLOCATE(p)
      
   ELSE
      PRINT*,'exiting now'
      STOP
   ENDIF
ENDIF

nullify(pt,rvalues,columns,rowIndex,div,p)
1000 FORMAT ('l2 of residual value is',e14.4)
1001 FORMAT ('l2 of residual value after is',e14.4)

END SUBROUTINE pardiso_solution





SUBROUTINE project(Info,cbc)
! subtract the gradient of the pressure field from the edge velocities
LOGICAL cbc
TYPE(Elli) :: Info
REAL (KIND=Precision), POINTER, DIMENSION (:,:,:) :: dummy,p,q,q1
INTEGER i,nx,ny,nz,j,k,BCI,lb,ub
INTEGER, DIMENSION (3,3) :: d
REAL (KIND=Precision) :: dx,dy,dz,x,y,z,xi,eta,zeta
d(1:3,1)=(/1,0,0/);d(1:3,2)=(/0,1,0/);d(1:3,3)=(/0,0,1/)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
dx=Info%dX(1);dy=Info%dX(2);dz=Info%dX(3)
p=>Info%p
IF(cbc) BCI=APPLY_BC(Info,AFTER_POISSON) ! apply bc
DO i=1,Info%nDim
   dummy=>adv(Info%G(i)%crsptr,p,(/-2,-2,-2/)+d(:,i)) ! gradient of pressure
   lb=0;ub=1
   IF(Info%bc((i-1)*2+1)%bc_type==NM) lb=1
   IF(Info%bc((i*2  ))%bc_type==NM) ub=0
!   lb=1;ub=0
   Info%qe(1+lb*d(1,i):Info%mX(1)+ub*d(1,i),&
        1+lb*d(2,i):Info%mX(2)+ub*d(2,i),1+lb*d(3,i):Info%mX(3)+ub*d(3,i),i)= &
        Info%qe(1+lb*d(1,i):Info%mX(1)+ub*d(1,i),&
        1+lb*d(2,i):Info%mX(2)+ub*d(2,i),1+lb*d(3,i):Info%mX(3)+ub*d(3,i),i)- &
        dummy(1+lb*d(1,i):Info%mX(1)+ub*d(1,i),&
        1+lb*d(2,i):Info%mX(2)+ub*d(2,i),1+lb*d(3,i):Info%mX(3)+ub*d(3,i))


! note how the velocities on the boundaries are not updated if the boundary
! is NM

DEALLOCATE(dummy);nullify(dummy)
END DO

END SUBROUTINE project


SUBROUTINE take_div(Info,cbc)

TYPE(Elli) :: Info
LOGICAL :: cbc
REAL (KIND=Precision), POINTER, DIMENSION (:,:,:) :: div,u,v,w,p
REAL (KIND=Precision), POINTER, DIMENSION (:,:) :: div2d,b,p0
!REAL (KIND=Precision), DIMENSION (SIZE(Info%bc(5)%p,1)-1,SIZE(Info%bc(5)%p,2)-1) :: p0

INTEGER ::i,j,k,nx,ny,nz,stage,l,ierr
INTEGER, DIMENSION (3,3) :: d
INTEGER, DIMENSION (2,6) :: bnx,bny,bnz
REAL (KIND=Precision) :: dx,dy,dz,betax,betay
REAL (KIND=Precision) :: x,y,z,xi,eta,zeta
REAL (KIND=Precision) :: fl(4)
! calculates the divergence and applies correction to divergence
! incorporating the boundary conditions of the problem
! NOTE THAT bc(5)%p has to include ghost layers even for NM conditions
d(1:3,1)=(/1,0,0/);d(1:3,2)=(/0,1,0/);d(1:3,3)=(/0,0,1/)


nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
dx=Info%dX(1);dy=Info%dX(2);dz=Info%dX(3)


   

IF(ASSOCIATED(Info%div)) THEN
   DEALLOCATE(INFO%div);nullify(Info%div)
ENDIF

ALLOCATE(Info%div(Info%mx(1),Info%mX(2),Info%mX(3)))
Info%div=zero
IF(cbc) THEN ! apply the NM conditions to the velocity
   DO l=1,2*Info%nDim
      IF(Info%bc(l)%bc_type==NM) THEN
         i=(l-1)/2+1
         u=>Info%qe(1:Info%mX(1)+d(1,i),1:Info%mX(2)+d(2,i),1:Info%mX(3)+d(3,i),i)
         SELECT CASE (l)
         CASE (1,2)
            
            u(1+l/2*nx,:,:)=u(1+l/2*nx,:,:)-Info%bc(l)%p(1:ny,1:nz)
!            Info%bc(l)%p(:,:)=zero ! so that later is not reused
            
         CASE (3,4)
            u(:,1+l/4*ny,:)=u(:,1+l/4*ny,:)-Info%bc(l)%p(1:nx,1:nz)
!            Info%bc(l)%p(:,:)=zero ! so that later is not reused
            
         CASE (5,6)
            u(:,:,1+l/6*nz)=u(:,:,1+l/6*nz)-Info%bc(l)%p(1:nx,1:ny)
!            Info%bc(l)%p(:,:)=zero ! so that later is not reused
            
         END SELECT
      ENDIF
   END DO
ENDIF
! calculates the divergence within the interior points
DO i=1,Info%nDim
   
   Info%div=Info%div+&
        (Info%qe(1+d(1,i):Info%mX(1)+d(1,i),1+d(2,i):Info%mX(2)+d(2,i),1+d(3,i):Info%mX(3)+d(3,i),i)-&
        Info%qe(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),i))/Info%dX(i)
!!$   div=>adv(Info%G(i+3)%crsptr,Info%qe(1:Info%mX(1)+d(1,i),1:Info%mX(2)+d(2,i),1:Info%mX(3)+d(3,i),i),-d(:,i))
!!$
!!$   Info%div=Info%div+div
!!%  DEALLOCATE(div);NULLIFY(div)
ENDDO

IF(.NOT.cbc) RETURN ! do not modify div, useful for diagnostic

! zero the interior points in the pressure
p=>Info%p(1:nx,1:ny,1:nz)
p=zero
!!$
!!$Info%p=zero
ierr=APPLY_BC(Info,BEFORE_POISSON)

! apply the boundary conditions


DO i=1,Info%ndim
   u=>adv(Info%G(i)%crsptr,Info%p,(/-2,-2,-2/)+d(:,i))
   ! zero the edge velocities at the boundary if boundary is NM
   SELECT CASE (i)
   CASE (1)
    IF(Info%bc((i-1)*2+1)%bc_type==NM) u(1,1:Info%mX(2),1:Info%mX(3))=zero
   IF(Info%bc((i*2  ))%bc_type==NM)&
        u(Info%mX(1)+1,1:Info%mX(2),1:Info%mX(3))=zero
   CASE (2)
    IF(Info%bc((i-1)*2+1)%bc_type==NM) u(1:Info%mX(1),1,1:Info%mX(3))=zero
   IF(Info%bc((i*2  ))%bc_type==NM)&
        u(1:Info%mX(1),Info%mX(2)+1,1:Info%mX(3))=zero

   CASE (3)
    IF(Info%bc((i-1)*2+1)%bc_type==NM) u(1:Info%mX(1),1:Info%mX(2),1)=zero
   IF(Info%bc((i*2  ))%bc_type==NM)&
        u(1:Info%mX(1),1:Info%mX(2),Info%mX(3)+1)=zero

   END SELECT
!!$   div=>adv(Info%G(i+3)%crsptr,u,-d(:,i))
!!$   Info%div=Info%div-div

   Info%div=Info%div-&
        (u(1+d(1,i):Info%mX(1)+d(1,i),1+d(2,i):Info%mX(2)+d(2,i),1+d(3,i):Info%mX(3)+d(3,i))-&
        u(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)))/Info%dX(i)


DEALLOCATE(U);
!!$DEALLOCATE(div)
ENDDO
END SUBROUTINE take_div


FUNCTION slice2d(Info,nxl,nxu,nyl,nyu,nzl,nzu)
! given pointer p, return a pointer to a 2D subset of p 
! (indicated by n_l=n_u) 

REAL (KIND=Precision), DIMENSION (:,:), POINTER :: slice2d

TYPE(Elli) :: Info

INTEGER, INTENT (IN) :: nxl,nxu,nyl,nyu,nzl,nzu



IF(nxu==nxl) THEN
   slice2d=>Info%p(nxu,    nyl:nyu,nzl:nzu) ! +1 because here indexes
                                                 ! start from one
ELSEIF(nyl==nyu) THEN
   slice2d=>Info%p(nxl:nxu,nyu    ,nzl:nzu)
ELSEIF(nzu==nzl) THEN
   slice2d=>Info%p(nxl:nxu,nyl:nyu,nzl    )
ELSE 
   PRINT*,'slice2d is confused, there are no singleton dimensions'
   PRINT*,'dimension of input is'
   PRINT*,nxl,nxu,nyl,nyu,nzl,nzu
   STOP
ENDIF



END FUNCTION slice2d



LOGICAL FUNCTION check(Info,argument)

TYPE (Elli), INTENT (IN) :: Info
INTEGER, INTENT (IN) :: argument

INTEGER :: nx,ny,nz,i,j,l,d(2)

nx=Info%mX(1);ny=Info%mX(2);nz=Info%mx(3)
check=.true.
SELECT CASE (argument)

CASE (CHECK_BC)
! check that allocation of bc is consistent

   DO i=1,Info%nDim
      SELECT CASE (i)
      CASE (1)
         d=(/Info%mx(2),Info%mX(3)/)
      CASE (2)
         d=(/Info%mX(1),Info%mX(3)/)
      CASE (3)
         d=(/Info%mX(1),Info%mX(2)/)
      CASE DEFAULT
         PRINT*,'If you get this message, you have some'
         PRINT*,'serious stack corruption problem...'
         STOP
      END SELECT
      DO j=1,2
         l=(i-1)*2+j
         SELECT CASE (l)
         CASE (5,6)
            check=(check.AND.ALL(SHAPE(Info%bc(l)%p).EQ.(d+2)))
            IF(.NOT.CHECK) THEN
               PRINT*,'BC ',l
               PRINT*,'is not consistent'
               PRINT*,'shape of array is'
               PRINT*, SHAPE(Info%bc(l)%p)
               PRINT*,'while it should be'
               PRINT*,d+2
               RETURN
            ENDIF
         CASE DEFAULT
            IF(Info%bc(l)%bc_type==NM) THEN        
               
               check=(check.AND.ALL(SHAPE(Info%bc(l)%p).EQ.(d+2)))
               IF(.NOT.CHECK) THEN
                  PRINT*,'BC ',l
                  PRINT*,'is not consistent'
                  PRINT*,'shape of array is'
                  PRINT*, SHAPE(Info%bc(l)%p)
                  PRINT*,'while it should be'
                  PRINT*,d+2
                  RETURN
               ENDIF
               
            ELSE
               check=(check.AND.ALL(SHAPE(Info%bc(l)%p).EQ.(d+2)))
               IF(.NOT.CHECK) THEN
                  PRINT*,'BC ',l
                  PRINT*,'is not consistent'
                  PRINT*,'shape of array is'
                  PRINT*, SHAPE(Info%bc(l)%p)
                  PRINT*,'while it should be'
                  PRINT*,d+2
                  RETURN
               ENDIF
               
            ENDIF
         END SELECT
      ENDDO
   ENDDO
   

                  
                  
                  
CASE DEFAULT
   PRINT*,'CASE ',argument
   PRINT*,'NOT IMPLEMENTED YET IN check'
   STOP
END SELECT
END FUNCTION check


INTEGER FUNCTION PARDISO_SOLVE(Info,action)
TYPE(Elli) :: Info
INTEGER, INTENT(IN) :: action
!
TYPE (CRS), POINTER :: A=>null()
SELECT CASE (action)

CASE (INIT)
CALL setparm(Info,PARDISO_SOLVER)
! extract symmetric part from LAPLACIAN
A=>EXTRACT_U_A(Info%G(CSR_LAPL)%crsptr)

CALL FREEUPCRS(Info%G(CSR_LAPL)%crsptr)
Info%G(CSR_LAPL)%crsptr=>A


CALL pardiso_init(Info)

CASE (SOLVE)

   CALL pardiso_solution(Info)

CASE (FINISH)

   CALL pardiso_mopup(Info)

   DEALLOCATE(Info%it%ipar_pd)

CASE DEFAULT 
PRINT*,'action ',action
PRINT*,'not implemented in PARDISO_SOLVE'
STOP
END SELECT
PARDISO_SOLVE=ErrOk
END FUNCTION PARDISO_SOLVE

 INTEGER FUNCTION DCG_SOLVE(Info,action)

TYPE(Elli) :: Info
INTEGER, INTENT(IN) :: action
!
TYPE (CRS), POINTER :: A=>NULL()
REAL (KIND=Precision),  DIMENSION (Info%mX(1),Info%mX(2),Info%mX(3),4),TARGET :: tmp
INTEGER, PARAMETER :: RCG_MULTIPLY = 1, RCG_STOPPING_TEST = 2,RCG_PRECONDITION = 3,RCG_OK=0
REAL (KIND=Precision), DIMENSION (:,:,:), POINTER :: div,p
REAL (KIND=Precision), DIMENSION (:), POINTER :: dpar
INTEGER, DIMENSION (:),POINTER  :: ipar
INTEGER :: error
CHARACTER, DIMENSION(3):: MATDES
INTEGER :: N,RCI_Request,nx,ny,nz

MATDES(1)='D'
MATDES(2)='L'
MATDES(3)='N'
nx=Info%mX(1);ny=Info%mx(2);nz=Info%mx(3)
   N=PRODUCT(Info%mX(1:3))
SELECT CASE (action)

CASE (INIT)
! extract symmetric part from LAPLACIAN
! this step is possibly redundant, but just in case...
   A=>EXTRACT_U_A(Info%G(CSR_LAPL)%crsptr)
   CALL FREEUPCRS(Info%G(CSR_LAPL)%crsptr)
   Info%G(CSR_LAPL)%crsptr=>A


   ALLOCATE(Info%it%fpar(128))
   ALLOCATE(Info%it%ipar_it(128))

   dpar=>Info%it%fpar
   ipar=>Info%it%ipar_it
   CALL dcg_init(N,Info%p(1:nx,1:ny,1:nz),Info%div,RCI_request,ipar,dpar,tmp)
IF(RCI_request/=RCG_OK) THEN
   PRINT*,'Failed initialization of CG'
   STOP
ENDIF

! set the parameters for the iteration here
! residual stopping test
! no user stopping test
! don't use preconditioner
! set tolerance

ipar(9:11)=(/1,0,1/);dpar(1)=1.d-18
! set max # of iteration to 1000
ipar(5)=5
CALL dcg_check(N,Info%p(1:nx,1:ny,1:nz),Info%div,RCI_request,ipar,dpar,tmp)
IF(RCI_request/=RCG_OK) THEN
   PRINT*,'Failed check of CG'
   PRINT*,'RCI_request = ',RCI_request
   STOP
ENDIF

CASE (SOLVE)

RCI_request=-999
DO WHILE(RCI_request/=RCG_OK)

   CALL dcg(N,Info%p(1:nx,1:ny,1:nz),Info%div,RCI_request,ipar,dpar,tmp)

!!$PRINT*,'Request is ',RCI_request
!!$PRINT*,'norm of p is ',SUM(Info%p(1:nx,1:ny,1:nz)**2)
   IF(RCI_request<0) THEN
      PRINT*,'dcg encoutered a problem and exited abnormally with code'
      PRINT*,RCI_request
      RCI_request=RCG_OK
   ENDIF

   SELECT CASE (RCI_request)
   CASE(RCG_MULTIPLY)
      !multiply A*tmp(:,1)=tmp(:,2)

      CALL mkl_dcsrsymv('U',N,Info%G(CSR_LAPL)%crsptr%rValues, &
           Info%G(CSR_LAPL)%crsptr%rowIndex,Info%G(CSR_LAPL)%crsptr%columns,tmp(1,1,1,1),tmp(1,1,1,2))

   CASE(RCG_STOPPING_TEST)
      PRINT*,'this is not required'
   CASE(RCG_PRECONDITION)

! to use the vertical derivative as preconditioner uncomment these lines

! save the original div and p

!!$      div=>Info%div
!!$      p=>Info%p
!!$      NULLIFY(Info%div,Info%p)
!!$      ALLOCATE(Info%div(Info%mX(1),Info%mX(2),Info%mX(3)))
!!$      ALLOCATE(Info%p(Info%mX(1),Info%mX(2),Info%mX(3)))
!!$
!!$      Info%div=RESHAPE(tmp(1:n,3),(/Info%mX(1),Info%mX(2),Info%mX(3)/))
!!$      Info%p(1:nx,1:ny,1:nz)=RESHAPE(tmp(1:n,4),(/Info%mX(1),Info%mX(2),Info%mX(3)/))
!!$      Info%div=>tmp(:,:,:,3);Info%p=>tmp(:,:,:,4)
!!$      error=POISSON_TRIDIAG(Info)
!!$      tmp(:,4:4)=RESHAPE(Info%p(1:nx,1:ny,1:nz),(/n,1/))
!!$
!!$
!!$      DEALLOCATE(Info%p,Info%div)
!!$      ! restore p and div
!!$      NULLIFY(Info%div,Info%p)
!!$      Info%div=>div
!!$      Info%p=>p
!!$
!!$      NULLIFY(div,p)

! to use the jacobi (i.e. divide by main diagonal) uncomment these

!!$      CALL MKL_DCSRSV('N', N, ONE, MATDES,Info%G(CSR_LAPL)%crsptr%rValues, &
!!$           Info%G(CSR_LAPL)%crsptr%columns,Info%G(CSR_LAPL)%crsptr%rowIndex,&
!!$           Info%G(CSR_LAPL)%crsptr%rowIndex(2),&
!!$              TMP(1,1,1,3), TMP(1,1,1,4))

! to do nothing uncomment these
      tmp(:,:,:,4)=tmp(:,:,:,3)

   END SELECT
END DO
PRINT*,'solved with ',ipar(4),' iterations'
CASE (FINISH)

   DEALLOCATE(Info%it%ipar_it,Info%it%fpar)
END SELECT
DCG_SOLVE=ErrOk
END FUNCTION DCG_SOLVE

RECURSIVE FUNCTION POISSON_SOLVERS(Info,action,solver,OP,LHS,RHS) RESULT(err)

TYPE(Elli) :: Info
INTEGER, INTENT (IN) :: solver,action
INTEGER :: err
INTEGER, INTENT (IN), OPTIONAL :: OP
REAL (KIND=Precision), DIMENSION (:), OPTIONAL :: LHS,RHS
! wrapper to Y.S. iterative solvers.
! the optional arguments and the recursive nature allow to use the
! direct solver as a preconditioner for the iterative solvers
! To do so, during initialize call with OP=position in Info%G
! where the preconditioner is stored
! and to solve the preconditioning system call with OP,LHS and RHS present
external  ilut


REAL (KIND=Precision), ALLOCATABLE, DIMENSION (:) ::w
INTEGER, DIMENSION (2*SIZE(Info%G(CSR_LAPL)%crsptr%rowIndex)-1) :: jw
INTEGER :: lfil,iwk,ierr=1,m,n
REAL (KIND=Precision) :: droptol
TYPE (CRS), POINTER :: A
TYPE (CRS), TARGET :: B
REAL (KIND=Precision), POINTER, DIMENSION (:) :: fpar
INTEGER, POINTER, DIMENSION (:) :: ipar
INTEGER (KIND=INT_PTR_KIND()), POINTER :: Addr

n=SIZE(Info%G(CSR_LAPL)%crsptr%rowIndex)-1

A=>Info%G(CSR_LAPL)%crsptr
SELECT CASE (action)

CASE (INIT)
   SELECT CASE (solver)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   CASE (DCG_SOLVER)
      err=DCG_SOLVE(Info,INIT)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   CASE (PARDISO_SOLVER)
      CALL setparm(Info,PARDISO_SOLVER)
      IF(.NOT.PRESENT(OP)) THEN
         ! extract symmetric part from LAPLACIAN
         ! the reason it is not done during generation
         ! is that ilut needs the full operator
         ! so in case ilut is to be used, then operator
         ! is not symmetrized when define, but needs to be done here
         A=>EXTRACT_U_A(Info%G(CSR_LAPL)%crsptr)
         CALL FREEUPCRS(INfo%G(CSR_LAPL)%crsptr)
         Info%G(CSR_LAPL)%crsptr=>A
         CALL pardiso_init(Info)
      ELSE
         CALL pardiso_init(Info,oper=OP)
      ENDIF
         
   CASE DEFAULT ! the iterative solvers have a common interface
      CALL setparm(Info,IT)
      IF(.NOT.PRESENT(OP)) THEN ! use ilut to generate preconditioner
         ALLOCATE(w(n+1))
      
         IF(ASSOCIATED(Info%it%G%rValues)) THEN ! clear up
            DEALLOCATE(Info%it%G%rValues); NULLIFY(Info%it%G%rValues)
         ENDIF
         IF(ASSOCIATED(Info%it%G%columns)) THEN ! clear up
            DEALLOCATE(Info%it%G%columns); NULLIFY(Info%it%G%columns)
         ENDIF
         IF(ASSOCIATED(Info%it%G%rowIndex)) THEN ! clear up
            DEALLOCATE(Info%it%G%rowIndex); NULLIFY(Info%it%G%rowIndex)
         ENDIF
         ! initial allocation, may change as needed
         iwk=SIZE(A%rValues) ! initial guess
         ALLOCATE(B%rValues(iwk))
         ALLOCATE(B%columns(iwk))
         ALLOCATE(B%rowIndex(n))
         
         
         droptol=Info%it%droptol
         lfil=Info%it%lfil
         DO WHILE(ierr/=0)
            call ilut (n,A%rValues,A%columns,A%rowIndex,lfil,droptol, &
                 B%rValues,B%columns,B%rowIndex,iwk,w,jw,ierr)
            
!!$            call milu0(n,A%rValues,A%columns,A%rowIndex, &
!!$                 B%rValues,B%columns,B%rowIndex,jw,ierr)
            SELECT CASE (ierr)
            CASE (1:)
               
               PRINT*,'Zero Pivot encountered in the generation of'
               PRINT*,'preconditioner at step, ',ierr
               PRINT*,'Stopping now'
               STOP
            CASE (-1) 
               PRINT*,'THe input matrix may be wrong or corrupted'
               PRINT*,'Stopping now'
               STOP
            CASE (-2,-3)
               DEALLOCATE(B%rvalues,B%columns)
               iwk=2*iwk
               ALLOCATE(B%rvalues(iwk),B%columns(iwk))
            CASE (-4)
               PRINT*,'Illegal value of ifil'
               PRINT*,'Fix and Restart'
               STOP
            CASE (-5) 
               PRINT*,'Zero row Encountered'
               PRINT*,'Stopping now'
               STOP
            END SELECT
         END DO
         
         Info%it%G=B
         DEALLOCATE(w);NULLIFY(A)
         ELSE
            err=POISSON_SOLVERS(Info,INIT,PARDISO_SOLVER,OP=OP)
            ! initialize the preconditioner
         ENDIF
         ! extract symmetric part from LAPLACIAN
         ! for rationale, see above
         A=>EXTRACT_U_A(Info%G(CSR_LAPL)%crsptr)
         CALL FREEUPCRS(INfo%G(CSR_LAPL)%crsptr)
         Info%G(CSR_LAPL)%crsptr=>A

      
      END SELECT
CASE (SOLVE)

   SELECT CASE (solver)
!!!!!!!!!!!!!!!!!!!!!!!!!!!
   CASE (DCG_SOLVER)
      err=DCG_SOLVE(Info,SOLVE)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
   CASE (PARDISO_SOLVER)
      IF(PRESENT(OP)) THEN
         CALL pardiso_solution(Info,oper=OP,rhs=RHS,lhs=LHS)
      ELSE
         CALL pardiso_solution(Info)
      ENDIF
      CASE DEFAULT ! the iterative solver
         fpar=>Info%it%fpar;ipar=>Info%it%ipar_it
         SELECT CASE (solver)
         CASE (IT_CG,IT_CGNR,IT_SD)
            iwk=5*n
            Info%it%ipar_it(4)=iwk
            ALLOCATE(w(iwk))





         CASE (IT_BCG)
            iwk=7*n
            Info%it%ipar_it(4)=iwk
            ALLOCATE(w(iwk))



         CASE (IT_DBCG)
            iwk=11*n
            Info%it%ipar_it(4)=iwk
            ALLOCATE(w(iwk))

         CASE (IT_BCGSTAB)
            iwk=8*n
            Info%it%ipar_it(4)=iwk
            ALLOCATE(w(iwk))

            
         CASE (IT_TFQMR)
            iwk=11*n
            Info%it%ipar_it(4)=iwk
            ALLOCATE(w(iwk))

         CASE (IT_FOM)
            m=ipar(5)
            iwk=(n+3)*(m+2) + (m+1)*m/2
            Info%it%ipar_it(4)=iwk
            ALLOCATE(w(iwk))

         CASE(IT_GMRES)
            m=ipar(5)
            iwk=(n+3)*(m+2) + (m+1)*m/2
            Info%it%ipar_it(4)=iwk

            ALLOCATE(w(iwk))

         CASE (IT_FGMRES)
            m=ipar(5)
            iwk=2*n*(m+1) + (m+1)*m/2 + 3*m + 2 
            Info%it%ipar_it(4)=iwk
            ALLOCATE(w(iwk))


         CASE (IT_DQGMRES)
            m=ipar(5)+1
            iwk=n + m * (2*n+4)
            Info%it%ipar_it(4)=iwk 
            ALLOCATE(w(iwk))

         END SELECT

         IF(.NOT.ASSOCIATED(Info%ptprcond)) THEN ! use ilut precond
            CALL runrc(Info,product(Info%mX(1:3)),&
                 Info%div,Info%p(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)), &
                 ipar,fpar,w,Info%pold(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)), &
                 A%rValues,A%columns,A%rowIndex, solver,&
                 au=Info%it%G%rValues,jau=Info%it%G%columns, &
                 ju=Info%it%G%rowIndex) 
         ELSE ! use pardiso to precondition
            CALL runrc(Info,product(Info%mX(1:3)),&
                 Info%div,Info%p(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)), &
                 ipar,fpar,w,Info%pold(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)), &
                 A%rValues,A%columns,A%rowIndex, &
                 solver)
            
         ENDIF
      DEALLOCATE(w)

      END SELECT
CASE (FINISH)
SELECT CASE (solver)
!!!!!!!!!!!!!!!!!!!!!!!
CASE (DCG_SOLVER)
   err=DCG_SOLVE(Info,FINISH)
!!!!!!!!!!!!!!!!!!!!!!!
CASE (PARDISO_SOLVER)
   IF(PRESENT(OP)) THEN ! used as a preconditioner
      CALL pardiso_mopup(Info,oper=OP)
   ELSE ! used as primary solver
      CALL pardiso_mopup(Info)

   ENDIF
      DEALLOCATE(Info%it%ipar_pd)

CASE DEFAULT
   IF(ASSOCIATED(Info%ptprcond)) THEN ! PARDISO preconditioner 
      err=POISSON_SOLVERS(Info,FINISH,PARDISO_SOLVER,OP=OP)
   ELSE
      DEALLOCATE(Info%it%G%rValues,Info%it%G%columns,Info%it%G%rowIndex)
   ENDIF
   DEALLOCATE(Info%it%fpar,Info%it%ipar_it)
   
END SELECT
END SELECT

err=ErrOk   
END FUNCTION POISSON_SOLVERS

subroutine runrc(Info,n,rhs,sol,ipar,fpar,wk,guess,a,ja,ia,solver,&
          au,jau,ju)
  TYPE(Elli) :: Info

  integer :: n,ipar(16),ia(n+1),err
  INTEGER, DIMENSION (:) :: ja
  INTEGER, DIMENSION (:), OPTIONAL :: ju,jau
  real (KIND=Precision) fpar(16),rhs(n),sol(n),guess(n)
  REAL (KIND=Precision), dimension (:) :: wk,a
  REAl (KIND=Precision), DIMENSION (:), OPTIONAL :: au
  INTEGER solver

  REAL (KIND=Precision), DIMENSION (:,:,:), POINTER :: div,p
  external mkl_dcsrsymv
EXTERNAL cg,bcg,dbcg,bcgstab,tfqmr,gmres,fgmres,dqgmres
EXTERNAL cgnr, fom
!!$c-----------------------------------------------------------------------
!!$c     the actual tester. It starts the iterative linear system solvers
!!$c     with a initial guess suppied by the user.
!!$c
!!$c     The structure {au, jau, ju} is assumed to have the output from
!!$c     the ILU* routines in ilut.f.
!!$c
!!$c-----------------------------------------------------------------------
!!$c     local variables
!!$c
      integer i, iou, its
      real*8 res, dnrm2
      save its,res
!!$c
!!$c     ipar(2) can be 0, 1, 2, please don't use 3
!!$c
      if (ipar(2).gt.2) then
         print *, 'I can not do both left and right preconditioning.'
         return
      endif
!!$c
!!$c     normal execution
!!$c
      its = 0
      res = 0.0D0
!c
      do i = 1, n
         sol(i) = guess(i)
      enddo
!c

      iou = 6
      ipar(1) = 0

 10   SELECT CASE (solver)
         CASE (IT_SD)
            call sd(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_CG)
            call cg(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_BCG)
            call bcg(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_DBCG)
            call dbcg(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_BCGSTAB)
            call bcgstab(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_TFQMR)
            CALL tfqmr(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_FOM)
            call fom(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_GMRES)

            call gmres(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_fgmres)
            call fgmres(n,rhs,sol,ipar,fpar,wk)
         CASE (IT_DQGMRES)

            call dqgmres(n,rhs,sol,ipar,fpar,wk)
         CASE DEFAULT
            PRINT*,'this solver ',solver
            PRINT*,'is not implemented'
            STOP
         END SELECT
!!$c
!!$c     output the residuals
!!$c
      if (ipar(7).ne.its) then
!         write (iou, *) its, real(res)
         its = ipar(7)
      endif
      res = fpar(5)

      if (ipar(1).eq.1) then

!         call amux(n, wk(ipar(8)), wk(ipar(9)), a, ja, ia)
         CALL mkl_dcsrsymv('U',N,a,ia,ja,wk(ipar(8)),wk(ipar(9)))
         goto 10
      else if (ipar(1).eq.2) then
         PRINT*,'stop'
         STOP
         call atmux(n, wk(ipar(8)), wk(ipar(9)), a, ja, ia)

         goto 10
      else if (ipar(1).eq.3 .or. ipar(1).eq.5) then


         IF(.NOT.PRESENT(au)) THEN

            err=POISSON_SOLVERS(Info,SOLVE,PARDISO_SOLVER,OP=CSR_PRECOND,rhs=wk(ipar(8):ipar(8)+N-1),lhs=wk(ipar(9):ipar(9)+N-1))

! the following does not seem to produce good results...

!!$         div=>Info%div;p=>Info%p
!!$         NULLIFY(Info%div,Info%p)
!!$         ALLOCATE(Info%div(Info%mX(1),Info%mX(2),Info%mX(3)))
!!$         ALLOCATE(Info%p(Info%mX(1),Info%mX(2),Info%mX(3)))
!!$         Info%div=RESHAPE(wk(ipar(8):ipar(8)+N-1),(/Info%mX(1),Info%mX(2),Info%mX(3)/))
!!$         Info%div(:,:,1)=Info%div(:,:,1)-SPREAD(SUM(Info%div(:,:,1),2)/Info%mX(2),DIM=2,NCOPIES=Info%mX(2))
!!$         err=POISSON_TRIDIAG(Info)
!!$
!!$         wk(ipar(9):ipar(9)+N-1)=RESHAPE(Info%p,(/N/))
!!$         DEALLOCATE(Info%p,Info%div)
!!$         Info%p=>p;Info%div=>div
!!$         NULLIFY(p,div)

         ELSE
         call lusol(n,wk(ipar(8)),wk(ipar(9)),au,jau,ju)
      ENDIF

         goto 10
      else if (ipar(1).eq.4 .or. ipar(1).eq.6) then

         call lutsol(n,wk(ipar(8)),wk(ipar(9)),au,jau,ju)
         goto 10
      else if (ipar(1).le.0) then
         if (ipar(1).eq.0) then
            print *, 'Iterative sovler has satisfied convergence test.'
         else if (ipar(1).eq.-1) then
            print *, 'Iterative solver has iterated too many times.'
         else if (ipar(1).eq.-2) then
            print *, 'Iterative solver was not given enough work space.'
            print *, 'The work space should at least have ', ipar(4),&
                ' elements.'
         else if (ipar(1).eq.-3) then
            print *, 'Iterative solver is facing a break-down.'
         else
            print *, 'Iterative solver terminated. code =', ipar(1)
         endif
      endif
!c     time = dtime(dt)
      write (iou, *) ipar(7), real(fpar(6))
      write (iou, *) '# return code =', ipar(1),&
          '	convergence rate =', fpar(7)
!!$c     write (iou, *) '# total execution time (sec)', time
!!$c
!!$c     check the error
!!$c
!!$c$$$      call amux(n,sol,wk,a,ja,ia)
!!$c$$$      do i = 1, n
!!$c$$$         wk(n+i) = sol(i) -1.0D0
!!$c$$$         WRITE(30,*),i,sol(i),wk(i),rhs(i)
!!$c$$$         wk(i) = wk(i) - rhs(i)
!!$c$$$
!!$c$$$      enddo
!!$c$$$      write (iou, *) '# the actual residual norm is', dnrm2(n,wk,1)
!!$c$$$      write (iou, *) '# the error norm is', dnrm2(n,wk(1+n),1)
!!$c
      if (iou.ne.6) close(iou)
      return
    end subroutine runrc
!!$c-----end-of-runrc
!!$c-----------------------------------------------------------------------



INTEGER FUNCTION FREEUP_OPERATORS(Info)

TYPE (Elli) :: Info
INTEGER i

! deallocates memory associated to matrix operators

DO i=LBOUND(Info%G,1),UBOUND(Info%G,1)
   IF(ASSOCIATED(Info%G(i)%crsptr)) &
        CALL FREEUPCRS(Info%G(i)%crsptr)
END DO
DO i=1,2*Info%nDim
   DEALLOCATE(info%bc(i)%p)
END DO

DEALLOCATE(Info%G);NULLIFY(Info%G)
DEALLOCATE(Info%bc);NULLIFY(Info%bc)
DEALLOCATE(Info%div);NULLIFY(info%div)
DEALLOCATE(Info%p);NULLIFY(Info%p)
DEALLOCATE(Info%pold);NULLIFY(Info%pold)
DEALLOCATE(Info%qe);nullify(Info%qe)
IF(ASSOCIATED(Info%diag)) DEALLOCATE(Info%diag);NULLIFY(Info%diag)

IF(ASSOCIATED(Info%it%fpar)) DEALLOCATE(Info%it%fpar)
IF(ASSOCIATED(Info%it%perm)) DEALLOCATE(Info%it%perm)
IF(ASSOCIATED(Info%it%ipar_pd)) DEALLOCATE(Info%it%ipar_pd)
IF(ASSOCIATED(Info%it%ipar_it)) DEALLOCATE(Info%it%ipar_it)
IF(ASSOCIATED(Info%it%G%rValues)) DEALLOCATE(Info%it%G%rValues)
IF(ASSOCIATED(Info%it%G%columns)) DEALLOCATE(Info%it%G%columns)
IF(ASSOCIATED(Info%it%G%rowIndex)) DEALLOCATE(Info%it%G%rowIndex)


FREEUP_OPERATORS=ErrOk
END FUNCTION FREEUP_OPERATORS


INTEGER FUNCTION APPLY_BC(Info,when)
! generates boundary conditions for the pressure

TYPE (Elli) :: Info
INTEGER, INTENT (IN) :: when
INTEGER, DIMENSION (3,3) :: d
INTEGER, DIMENSION (2,6) :: bnx,bny,bnz
REAL (KIND=Precision) :: dx,dy,dz,betax,betay
REAL (KIND=Precision) :: x,y,z,xi,eta,zeta
REAL (KIND=Precision) :: fl(4)
REAL (KIND=Precision), POINTER, DIMENSION (:,:) :: b,p0
REAL (KIND=Precision), DIMENSION (:,:,:), POINTER :: dummy
REAl (KIND=Precision), DIMENSION (:), ALLOCATABLE :: dl,pl
INTEGER :: i,j,k,l,nx,ny,nz,kc,jc,ic
d(1:3,1)=(/1,0,0/);d(1:3,2)=(/0,1,0/);d(1:3,3)=(/0,0,1/)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)


! fill in the boundaries arrays

IF(Info%bc(1)%bc_type==NM) THEN
   bnx(:,1)=0
   bny(:,1)=(/1,ny/)
   bnz(:,1)=(/1,nz/)
ELSE
   bnx(:,1)=0
   bny(:,1)=(/0,ny+1/)
   bnz(:,1)=(/0,nz+1/)
ENDIF

IF(Info%bc(2)%bc_type==NM) THEN
   bnx(:,2)=nx+1
   bny(:,2)=(/1,ny/)
   bnz(:,2)=(/1,nz/)
ELSE
   bnx(:,2)=nx+1
   bny(:,2)=(/0,ny+1/)
   bnz(:,2)=(/0,nz+1/)
ENDIF

IF(Info%bc(3)%bc_type==NM) THEN
   bnx(:,3)=(/1,nx/)
   bny(:,3)=0
   bnz(:,3)=(/1,nz/)
ELSE
   bnx(:,3)=(/0,nx+1/)
   bny(:,3)=0
   bnz(:,3)=(/0,nz+1/)

ENDIF

IF(Info%bc(4)%bc_type==NM) THEN
   bnx(:,4)=(/1,nx/)
   bny(:,4)=ny+1
   bnz(:,4)=(/1,nz/)
ELSE
   bnx(:,4)=(/0,nx+1/)
   bny(:,4)=ny+1
   bnz(:,4)=(/0,nz+1/)

ENDIF
IF(Info%nDim>2) THEN ! third dimension
   IF(Info%bc(5)%bc_type==NM) THEN
      bnx(:,5)=(/1,nx/) 
      bny(:,5)=(/1,ny/)
      bnz(:,5)=0
   ELSE
      bnx(:,5)=(/0,nx+1/) ! 
      bny(:,5)=(/0,ny+1/)
      bnz(:,5)=0
   ENDIF
   
   IF(Info%bc(6)%bc_type==NM) THEN
      bnx(:,6)=(/1,nx/) 
      bny(:,6)=(/1,ny/)
      bnz(:,6)=nz+1
   ELSE
      bnx(:,6)=(/0,nx+1/) !
      bny(:,6)=(/0,ny+1/)
      bnz(:,6)=nz+1
   ENDIF
ENDIF
! 1. calculates DR extrapolations based on BC

SELECT CASE (when)

   CASE (BEFORE_POISSON) 
      ! in this case only the nonhomogeneous values of BC
      ! are filled

      DO l=1,2*Info%ndim
         p0=>slice2d(Info,bnx(1,l),bnx(2,l),bny(1,l),bny(2,l),bnz(1,l),bnz(2,l))

         IF(Info%bc(l)%bc_type==DR) THEN
              p0=Info%bc(l)%p
         ELSE
            SELECT CASE (l)
            CASE(1,2)
               i=1+l/2*nx
               DO k=1,nz
                  DO j=1,ny

                     p0(j,k)=(-1)**l*Info%dX((l-1)/2+1)*Info%bc(l)%p(j,k)/g(i,j,k,1,1)

                  ENDDO
               ENDDO
            CASE (3,4)
               j=1+l/4*ny
               DO k=1,nz
                  DO i=1,nx
                     
                     p0(i,k)=(-1)**l*Info%dX((l-1)/2+1)*Info%bc(l)%p(i,k)/g(i,j,k,2,2)
                  ENDDO
               ENDDO
            CASE(5,6)
               
               k=1+l/6*nz
               DO j=1,ny
                  DO i=1,nx
                     
                     p0(i,j)=(-1)**l*Info%dX((l-1)/2+1)*&
                          Info%bc(l)%p(i,j)/g(i,j,k,3,3)
                  ENDDO
               ENDDO
            END SELECT


         ENDIF

       END DO


          
   CASE (AFTER_POISSON)

      dummy=>adv(Info%G(EXTR_BC)%crsptr,Info%p(1:nx,1:ny,1:nz),(/2,2,2/))

      Info%p(0:nx+1,0:ny+1,0:nz+1)=Info%p(0:nx+1,0:ny+1,0:nz+1)+dummy(1:nx+2,1:ny+2,1:nz+2)
      
      DEALLOCATE(dummy)
      
      ! now add the zeroth_order extrapolation 

   CASE DEFAULT
         
      PRINT*,'APPLY_BC has been called with an unknown option'
      PRINT*,'terminating now'
      STOP
         
      
         
   END SELECT
   APPLY_BC=ErrOk
 END FUNCTION APPLY_BC

INTEGER FUNCTION CREATE_OPERATORS(Info)

TYPE (Elli) :: Info
INTEGER :: nx,ny,nz,i,n,st,d(3,3),ierr,j,k,l
TYPE (CRS), POINTER :: crsptr=>NULL()
TYPE (CRS), POINTER :: B=>NULL(),A=>NULL(),C=>NULL(),DC=>NULL()
REAL(KIND=Precision), POINTER, DIMENSION (:,:,:) :: dummy

REAL(KIND=Precision), DIMENSION (:), POINTER :: rv
INTEGER, DIMENSION (:), POINTER :: ja
d(1:3,1)=(/1,0,0/);d(1:3,2)=(/0,1,0/);d(1:3,3)=(/0,0,1/)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
! used to create and allocate the operators used during the Projection
! step. 
! last allocate the space for the arrays

ALLOCATE(Info%p(0:nx+1,0:ny+1,0:nz+1))
ALLOCATE(Info%pold(1:nx,1:ny,1:nz))
ALLOCATE(Info%div(nx,ny,nz))
ALLOCATE(Info%qe(nx+1,ny+1,nz+1,Info%ndim))
Info%p=zero;Info%pold=zero;Info%div=zero;Info%qe=zero
CREATE_OPERATORS=ErrOk
! 0. Allocate the structures to hold the operators and the boundary conditions

ALLOCATE(Info%G(0:8))
DO i=0,8
   Info%G(i)%crsptr=>NULL() ! to avoid issues later
ENDDO
ALLOCATE(Info%bc(6))
! 1. create the laplacian

! 1.1 set the type of boundary conditions
DO i=1,2*Info%nDim
   IF(Info%mthbc(i)<InternalBC) THEN
      Info%bc(i)%bc_type=NM
   ELSE
      Info%bc(i)%bc_type=DR
   ENDIF
ENDDO

i=1

   ALLOCATE(Info%bc(i)%p(0:ny+1,0:nz+1))

i=i+1

   ALLOCATE(Info%bc(i)%p(0:ny+1,0:nz+1))

i=i+1
   ALLOCATE(Info%bc(i)%p(0:nx+1,0:nz+1))

i=i+1

   ALLOCATE(Info%bc(i)%p(0:nx+1,0:nz+1))
IF(Info%nDim>2) THEN
i=i+1

   ALLOCATE(Info%bc(i)%p(0:nx+1,0:ny+1))

i=i+1
   ALLOCATE(Info%bc(i)%p(0:nx+1,0:ny+1))
ENDIF

! operator to extrapolate BC (homogeneous case)
Info%G(EXTR_BC)%crsptr=>EXTRAPOLATE_BC(Info)

! divergence
DO i=4,4+Info%nDim-1
   ALLOCATE(dummy(nx+d(1,i-3),ny+d(2,i-3),nz+d(3,i-3)))
   Info%G(i)%crsptr=>operation('ge2cs',dummy,i-3,delta=Info%dX(i-3))
   DEALLOCATE(dummy)
ENDDO

DO i=1,Info%nDim ! create the gradient operators with DR conditions



Info%G(i)%crsptr=>GRAD(Info,i)



ENDDO

! generate laplacian
! store matrix for BC
   ALLOCATE(dummy(nx,ny,nz))
   B=>operation('dress',dummy,0)
   DEALLOCATE(dummy)
   DC=>APB(B,Info%G(EXTR_BC)%crsptr);CALL FREEUPCRS(B)
! x direction
   B=>ADB(Info%G(1)%crsptr,DC) ! apply bc

   ! zero gradient at neumann x-boundaries
   IF(Info%bc(1)%bc_type==NM) THEN
      FORALL(k=1:nz,j=1:ny)
         B%rValues(B%rowIndex((k-1)*ny*(nx+1)+(j-1)*(nx+1)+1):&
              B%rowIndex((k-1)*ny*(nx+1)+(j-1)*(nx+1)+1+1)-1)=zero
      END FORALL
   END IF
   
   IF(Info%bc(2)%bc_type==NM) THEN
      FORALL(k=1:nz,j=1:ny)
         
         B%rValues(B%rowIndex((k-1)*ny*(nx+1)+(j-1)*(nx+1)+nx+1):&
              B%rowIndex((k-1)*ny*(nx+1)+(j-1)*(nx+1)+nx+1+1)-1)=zero
      END FORALL

   END IF
A=>ADB(Info%G(4)%crsptr,B);CALL FREEUPCRS(B) !d1*g1j*dj

C=>ADB(Info%G(2)%crsptr,DC) ! apply bc

! zero gradient along y NM boundaries
IF(Info%bc(3)%bc_type==NM) THEN

   FORALL(k=1:nz, i=1:nx)
      
      C%rValues(C%rowIndex((k-1)*(ny+1)*nx+i):&
           C%rowIndex((k-1)*(ny+1)*nx+i+1)-1)=zero
      
   END FORALL
END IF

IF(Info%bc(4)%bc_type==NM) THEN
   FORALL (k=1:nz, i=1:nx)
      C%rValues(C%rowIndex((k-1)*(ny+1)*nx+ny*nx+i):&
           C%rowIndex((k-1)*(ny+1)*nx+ny*nx+i+1)-1)=zero
   END FORALL
   
END IF



B=>ADB(Info%G(5)%crsptr,C);CALL FREEUPCRS(C) ! d2*g2j*dj

IF(Info%nDim==2) THEN

   Info%G(CSR_LAPL)%crsptr=>APB(A,B); !d1*g1j*dj+d2*g2j*dj
   CALL FREEUPCRS(A);CALL FREEUPCRS(B)

ELSE ! add third direction
   C=>APB(A,B)!d1*g1j*dj+d2*g2j*dj
   CALL FREEUPCRS(A);CALL FREEUPCRS(B)

   A=>ADB(Info%G(3)%crsptr,DC) ! add bc

! zero along z NM planes
   IF(Info%bc(5)%bc_type==NM) THEN
      k=1
      FORALL(j=1:ny,i=1:nx)
         A%rValues(A%rowIndex((j-1)*nx+i):&
              A%rowIndex((j-1)*nx+i+1)-1)=zero
      END FORALL
   ENDIF

   IF(Info%bc(6)%bc_type==NM) THEN
      k=nz+1
      FORALL(j=1:ny,i=1:nx)
         A%rValues(A%rowIndex(nz*nx*ny+(j-1)*nx+i):&
              A%rowIndex(nz*nx*ny+(j-1)*nx+i+1)-1)=zero
      END FORALL
   ENDIF
   
 

   B=>ADB(Info%G(6)%crsptr,A);CALL FREEUPCRS(A) ! d3*g3j*dj
   
   Info%G(CSR_LAPL)%crsptr=>APB(C,B)!d1*g1j*dj+d2*g2j*dj+d3*g3j*dj

   CALL FREEUPCRS(C);CALL FREEUPCRS(B)
   

ENDIF
CALL FREEUPCRS(DC)
! filter small entries from laplacian
Info%G(CSR_LAPL)%crsptr=>FILTER_Z_P(Info%G(CSR_LAPL)%crsptr,1.d-14)



   

! create the operator for the preconditioning
! strategy is to zero all the diagonals in CSR_LAPL except the one corresponding
! to second derivative in z 

! get upper triangular

A=>EXTRACT_U_P(Info%G(CSR_LAPL)%crsptr,keep=.TRUE.)

n=SIZE(A%rowIndex)-1 ! # of rows
! remove diagonals

IF(Info%nDim==3) THEN
   st=nx*ny
ELSE
   st=nx
ENDIF
!!$FORALL (i=1:n)
DO i=1,n
   ja=>A%columns(A%rowIndex(i):A%rowIndex(i+1)-1)
   rv=>A%rValues(A%rowIndex(i):A%rowIndex(i+1)-1)

   WHERE((ja/=i).AND.(ja/=i+st))
      rv=zero ! zeros the off diagonal terms, except k+1 terms
   END WHERE
ENDDO
!!$END FORALL

! remove zero entries
Info%G(CSR_PRECOND)%crsptr=>FILTER_Z_P(A,1.d-14)


      



CALL FREEUPCRS(A);


END FUNCTION CREATE_OPERATORS

INTEGER FUNCTION CREATE_OPERATORS_THIN(Info)
! same functionality as create_operators, but
! assume thin layer approximation to reduce dimensionality
! of problem. See notes...
TYPE (Elli) :: Info
INTEGER :: nx,ny,nz,i,n,st,d(3,3),ierr,j,k,l,norig
TYPE (CRS), POINTER :: crsptr=>NULL()
TYPE (CRS), POINTER :: B=>NULL(),A=>NULL(),C=>NULL(),DC=>NULL()
REAL(KIND=Precision), POINTER, DIMENSION (:,:,:) :: dummy

REAL(KIND=Precision), DIMENSION (:), POINTER :: rv
INTEGER, DIMENSION (:), POINTER :: ja
!!$IF(Info%nDim<3) THEN
!!$   PRINT*,'2D problem not implemented yet for thin'
!!$   PRINT*,'case ...'
!!$   STOP
!!$ENDIF
d(1:3,1)=(/1,0,0/);d(1:3,2)=(/0,1,0/);d(1:3,3)=(/0,0,1/)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
! used to create and allocate the operators used during the Projection
! step. 
! last allocate the space for the arrays

ALLOCATE(Info%p(0:nx+1,0:ny+1,0:nz+1))
ALLOCATE(Info%pold(1:nx,1:ny,1:nz))
ALLOCATE(Info%div(nx,ny,nz))
ALLOCATE(Info%qe(nx+1,ny+1,nz+1,Info%ndim))

Info%p=zero;Info%pold=zero;Info%div=zero;Info%qe=zero

CREATE_OPERATORS_THIN=ErrOk
! 0. Allocate the structures to hold the operators and the boundary conditions

ALLOCATE(Info%G(0:10))
! scheme 0 -> extr_BC2d
! 1-3 -> grad (3D)
! 4-6 -> div (3d)
! 7 -> laplacian (2d)
! 8 -> d3g33d3 (3)
! 9 -> d3g31d1+d3g32d2
! 10 -> g31d1+g32d2
 
DO i=0,10
   Info%G(i)%crsptr=>NULL() ! to avoid issues later
ENDDO
ALLOCATE(Info%bc(6))
! 1. create the laplacian

! 1.1 set the type of boundary conditions
DO i=1,2*Info%nDim
   IF(Info%mthbc(i)<InternalBC) THEN
      Info%bc(i)%bc_type=NM
   ELSE
      Info%bc(i)%bc_type=DR
   ENDIF
ENDDO

i=1

   ALLOCATE(Info%bc(i)%p(0:ny+1,0:nz+1))

i=i+1

   ALLOCATE(Info%bc(i)%p(0:ny+1,0:nz+1))

i=i+1
   ALLOCATE(Info%bc(i)%p(0:nx+1,0:nz+1))

i=i+1

   ALLOCATE(Info%bc(i)%p(0:nx+1,0:nz+1))
IF(Info%nDim>2) THEN
i=i+1

   ALLOCATE(Info%bc(i)%p(0:nx+1,0:ny+1))

i=i+1
   ALLOCATE(Info%bc(i)%p(0:nx+1,0:ny+1))
ENDIF

! operator to extrapolate BC (homogeneous case)
Info%G(EXTR_BC3d)%crsptr=>EXTRAPOLATE_BC(Info)
! reduce dimension and setup 2D problem
Norig=Info%mX(Info%ndim)
Info%mX(Info%ndim)=1;Info%ndim=Info%ndim-1;Info%G(EXTR_BC2d)%crsptr=>EXTRAPOLATE_BC(Info)
nz=1
IF(Info%ndim==1) ny=1
! divergence (2D)
DO i=4,4+Info%nDim-1
   ALLOCATE(dummy(nx+d(1,i-3),ny+d(2,i-3),nz+d(3,i-3)))
   Info%G(i)%crsptr=>operation('ge2cs',dummy,i-3,delta=Info%dX(i-3))
   DEALLOCATE(dummy)
ENDDO
! gradient (2D)
DO i=1,Info%nDim ! create the gradient operators with DR conditions
Info%G(i)%crsptr=>GRAD(Info,i)
ENDDO

! generate laplacian 2d


! store matrix for BC
   ALLOCATE(dummy(nx,ny,nz))
   B=>operation('dress',dummy,0)
   DEALLOCATE(dummy)
   DC=>APB(B,Info%G(EXTR_BC2d)%crsptr);CALL FREEUPCRS(B)
! x direction
   B=>ADB(Info%G(1)%crsptr,DC) ! apply bc

   ! zero gradient at neumann x-boundaries
   IF(Info%bc(1)%bc_type==NM) THEN
      FORALL(k=1:nz,j=1:ny)
         B%rValues(B%rowIndex((k-1)*ny*(nx+1)+(j-1)*(nx+1)+1):&
              B%rowIndex((k-1)*ny*(nx+1)+(j-1)*(nx+1)+1+1)-1)=zero
      END FORALL
   END IF
   
   IF(Info%bc(2)%bc_type==NM) THEN
      FORALL(k=1:nz,j=1:ny)
         
         B%rValues(B%rowIndex((k-1)*ny*(nx+1)+(j-1)*(nx+1)+nx+1):&
              B%rowIndex((k-1)*ny*(nx+1)+(j-1)*(nx+1)+nx+1+1)-1)=zero
      END FORALL

   END IF
IF(Info%ndim==1) THEN

   Info%G(CSR_LAPL)%crsptr=>ADB(Info%G(4)%crsptr,B)
   CALL FREEUPCRS(B)
ELSE

   A=>ADB(Info%G(4)%crsptr,B);CALL FREEUPCRS(B) !d1*g1j*dj
   


   C=>ADB(Info%G(2)%crsptr,DC) ! apply bc

   ! zero gradient along y NM boundaries
   IF(Info%bc(3)%bc_type==NM) THEN
      
      FORALL(k=1:nz, i=1:nx)
         
         C%rValues(C%rowIndex((k-1)*(ny+1)*nx+i):&
              C%rowIndex((k-1)*(ny+1)*nx+i+1)-1)=zero
         
      END FORALL
   END IF
   
   IF(Info%bc(4)%bc_type==NM) THEN
      FORALL (k=1:nz, i=1:nx)
         C%rValues(C%rowIndex((k-1)*(ny+1)*nx+ny*nx+i):&
              C%rowIndex((k-1)*(ny+1)*nx+ny*nx+i+1)-1)=zero
      END FORALL
      
   END IF
   


   B=>ADB(Info%G(5)%crsptr,C);CALL FREEUPCRS(C) ! d2*g2j*dj



   Info%G(CSR_LAPL)%crsptr=>APB(A,B); !d1*g1j*dj+d2*g2j*dj
   CALL FREEUPCRS(A);CALL FREEUPCRS(B)
END IF
CALL FREEUPCRS(DC)

! filter small entries from laplacian
Info%G(CSR_LAPL)%crsptr=>FILTER_Z_P(Info%G(CSR_LAPL)%crsptr,1.d-14)




   

! create the operator for the vertical integration
Info%ndim=Info%ndim+1 ! restore dimensionality to original
Info%mX(Info%nDim)=nOrig
ny=Info%mX(2);nz=Info%mX(3)
! in diag we store |g|g^33 to be used by tridiag solver
IF(Info%nDim==2) THEN
   ALLOCATE(Info%diag(1:nx,1:ny+1,1:nz))
   FORALL (i=1:nx,j=1:ny+1)
      Info%diag(i,j,1)=g(i,j,1,2,2)
   END FORALL
ELSE
   ALLOCATE(Info%diag(1:nx,1:ny,1:nz+1))
   FORALL (i=1:nx,j=1:ny,k=1:nz+1)
      Info%diag(i,j,k)=g(i,j,k,3,3)
   END FORALL
END IF

!!$PRINT*,nx,ny,nz
!!$PRINT*,Info%ndim
!!$PRINT*,Info%mX
DO i=1,6
   CALL FREEUPCRS(Info%g(i)%crsptr)
ENDDO
! divergence (3D)
DO i=4,4+Info%nDim-1
   ALLOCATE(dummy(nx+d(1,i-3),ny+d(2,i-3),nz+d(3,i-3)))
   Info%G(i)%crsptr=>operation('ge2cs',dummy,i-3,delta=Info%dX(i-3))
   DEALLOCATE(dummy)
ENDDO
! gradient (3D)
DO i=1,Info%nDim ! create the gradient operators with DR conditions
Info%G(i)%crsptr=>GRAD(Info,i)
ENDDO


! create matrix for vertical problem
! store matrix for BC
!!$   ALLOCATE(dummy(nx,ny,nz))
!!$   B=>operation('dress',dummy,0)
!!$   DEALLOCATE(dummy)
!!$   DC=>APB(B,Info%G(EXTR_BC3d)%crsptr);CALL FREEUPCRS(B)
!!$   IF(Info%ndim==3) THEN
!!$! z direction
!!$      B=>gijdj(Info,3,3,'b')
!!$      A=>ADB(B,DC) ! apply bc
!!$      CALL FREEUPCRS(B);
!!$      IF(Info%bc(5)%bc_type==NM) THEN
!!$         k=1
!!$         FORALL(j=1:ny,i=1:nx)
!!$            A%rValues(A%rowIndex((j-1)*nx+i):&
!!$                 A%rowIndex((j-1)*nx+i+1)-1)=zero
!!$         END FORALL
!!$      ENDIF
!!$
!!$      IF(Info%bc(6)%bc_type==NM) THEN
!!$         k=nz+1
!!$         FORALL(j=1:ny,i=1:nx)
!!$            A%rValues(A%rowIndex(nz*nx*ny+(j-1)*nx+i):&
!!$                 A%rowIndex(nz*nx*ny+(j-1)*nx+i+1)-1)=zero
!!$         END FORALL
!!$      ENDIF
!!$      
!!$      B=>ADB(Info%G(6)%crsptr,A);CALL FREEUPCRS(A) ! d3*g33*d3
!!$      
!!$      Info%G(CSR_PRECOND)%crsptr=>FILTER_Z_P(B,1.d-14)
!!$      CALL FREEUPCRS(B)
!!$   ELSE
!!$      B=>gijdj(Info,2,2,'b')
!!$      C=>ADB(B,DC); CALL FREEUPCRS(B)
!!$      
!!$      IF(Info%bc(3)%bc_type==NM) THEN
!!$         
!!$         FORALL(k=1:nz, i=1:nx)
!!$            
!!$            C%rValues(C%rowIndex((k-1)*(ny+1)*nx+i):&
!!$                 C%rowIndex((k-1)*(ny+1)*nx+i+1)-1)=zero
!!$            
!!$         END FORALL
!!$      END IF
!!$      
!!$      IF(Info%bc(4)%bc_type==NM) THEN
!!$         FORALL (k=1:nz, i=1:nx)
!!$            C%rValues(C%rowIndex((k-1)*(ny+1)*nx+ny*nx+i):&
!!$                 C%rowIndex((k-1)*(ny+1)*nx+ny*nx+i+1)-1)=zero
!!$         END FORALL
!!$         
!!$      END IF
!!$      
!!$      
!!$      B=>ADB(Info%G(5)%crsptr,C);CALL FREEUPCRS(C) ! d3*g33*d3
!!$      Info%G(CSR_PRECOND)%crsptr=>FILTER_Z_P(B,1.d-14)
!!$      CALL FREEUPCRS(B)
!!$   END IF
   
   ! finally
!!$IF(Info%nDim==3) THEN
!!$   B=>gijdj(Info,3,1,'b')
!!$   A=>gijdj(Info,3,2,'b')
!!$   C=>APB(A,B); CALL FREEUPCRS(A); CALL FREEUPCRS(B)
!!$   A=>ADB(C,DC); CALL FREEUPCRS(C);CALL FREEUPCRS(DC)
!!$   B=>ADB(Info%G(6)%crsptr,A); CALL FREEUPCRS(A)
!!$   
!!$! remove zero entries
!!$   Info%G(9)%crsptr=>FILTER_Z_P(B,1.d-14)
!!$CALL FREEUPCRS(B);
!!$ELSE
!!$   B=>gijdj(Info,2,1,'b')
!!$   A=>ADB(B,DC);CALL FREEUPCRS(B);CALL FREEUPCRS(DC)
!!$   B=>ADB(Info%G(5)%crsptr,A);CALL FREEUPCRS(A)
!!$   Info%G(9)%crsptr=>FILTER_Z_P(B,1.d-14); CALL FREEUPCRS(B)
!!$ENDIF


      






END FUNCTION CREATE_OPERATORS_THIN

RECURSIVE FUNCTION POISSON_SOLVERS_THIN(Info,action,OP,LHS,RHS) RESULT(err)

TYPE(Elli) :: Info
INTEGER, INTENT (IN) :: action
INTEGER :: err
INTEGER, OPTIONAL :: OP
REAL (KIND=Precision), DIMENSION (:), OPTIONAL :: LHS,RHS
! wrapper to Y.S. iterative solvers.
! the optional arguments and the recursive nature allow to use the
! direct solver as a preconditioner for the iterative solvers
! To do so, during initialize call with OP=position in Info%G
! where the preconditioner is stored
! and to solve the preconditioning system call with OP,LHS and RHS present
external  ilut


REAL (KIND=Precision), ALLOCATABLE, DIMENSION (:) ::w
INTEGER, DIMENSION (2*SIZE(Info%G(CSR_LAPL)%crsptr%rowIndex)-1) :: jw
INTEGER :: lfil,iwk,ierr=1,m,n
REAL (KIND=Precision) :: droptol
TYPE (CRS), POINTER :: A
TYPE (CRS), TARGET :: B
REAL (KIND=Precision), POINTER, DIMENSION (:) :: fpar
INTEGER, POINTER, DIMENSION (:) :: ipar
INTEGER (KIND=INT_PTR_KIND()), POINTER :: Addr

n=SIZE(Info%G(CSR_LAPL)%crsptr%rowIndex)-1


SELECT CASE (action)

CASE (INIT)
      CALL setparm(Info,PARDISO_SOLVER)
      IF(.NOT.PRESENT(OP)) THEN
         ! extract symmetric part from LAPLACIAN
         ! the reason it is not done during generation
         ! is that ilut needs the full operator
         ! so in case ilut is to be used, then operator
         ! is not symmetrized when define, but needs to be done here
         A=>EXTRACT_U_A(Info%G(CSR_LAPL)%crsptr)
         CALL FREEUPCRS(INfo%G(CSR_LAPL)%crsptr)
         Info%G(CSR_LAPL)%crsptr=>A
         n=Info%mX(Info%ndim)
         Info%mX(Info%ndim)=1
         CALL pardiso_init(Info)
         Info%mX(Info%ndim)=n
         
         
      ELSE
         A=>EXTRACT_U_A(Info%G(OP)%crsptr)
         CALL FREEUPCRS(Info%G(OP)%crsptr)
         Info%G(OP)%crsptr=>A
         CALL pardiso_init(Info,oper=OP)
      ENDIF
         
CASE (SOLVE)

      IF(PRESENT(OP)) THEN
         CALL pardiso_solution(Info,oper=OP)
      ELSE
         n=Info%mX(Info%ndim)
         Info%mX(Info%ndim)=1
         
         CALL pardiso_solution(Info)
         Info%mX(Info%ndim)=n
      ENDIF

CASE (FINISH)


   IF(PRESENT(OP)) THEN
      CALL pardiso_mopup(Info,oper=OP)
   ELSE
      n=Info%mX(Info%ndim)
      Info%mX(Info%ndim)=1
      CALL pardiso_mopup(Info)
      Info%mX(Info%ndim)=n
   ENDIF
      DEALLOCATE(Info%it%ipar_pd)      
END SELECT

err=ErrOk   
END FUNCTION POISSON_SOLVERS_THIN

! ************************************************************

SUBROUTINE hodge_thin(Info,maxit,vs,ut,vt,wt)
!
TYPE (Elli) :: Info
REAL (KIND=Precision), INTENT (IN) :: ut(1:Info%mX(2),1:Info%mX(3),2)
REAL (KIND=Precision), INTENT (IN) :: vt(1:Info%mX(1),1:Info%mX(3),2)
REAL (KIND=Precision), OPTIONAL, INTENT (IN) :: wt(1:Info%mX(1),1:Info%mX(2),2)
INTEGER, INTENT (IN) :: maxit
LOGICAL, INTENT (IN) :: vs
!
REAL (KIND=Precision) :: pt(0:Info%mX(1)+1,0:Info%mX(2)+1,0:Info%mX(3)+1),&
     pacc(0:Info%mX(1)+1,0:Info%mX(2)+1,0:Info%mX(3)+1)
LOGICAL :: Iterate
INTEGER :: nx,ny,nz,ndim,N,i,error,j
REAL (KIND=Precision), DIMENSION (:,:,:), POINTER :: u,v,w,div,p,dummy,dummy2
REAL (KIND=Precision), DIMENSION (:,:,:), ALLOCATABLE :: p_temp
! sanity check

IF(.NOT.PRESENT(wt).AND.Info%ndim==3) THEN
   PRINT*,'wt must be input in 3D'
   STOP
END IF
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3);ndim=Info%ndim;N=Info%mX(ndim)

!!$pt=zero

pacc=zero
u=>Info%qe(1:nx+1,1:ny,1:nz,1)
v=>Info%qe(1:nx,1:ny+1,1:nz,2)
w=>Info%qe(1:nx,1:ny,1:nz+1,3)
div=>Info%div
p=>Info%p
!!$pt=p
pt=zero
Iterate=.true.;i=0

DO WHILE(Iterate.AND.i<maxit)

! *********************** vertical component *******************
! set boundary conditions
! left and right
Info%bc(1)%p(1:ny,1:nz)=u(1,:,:)-ut(:,:,1)-&
     SPREAD(SUM(u(1,:,:)-ut(:,:,1),ndim-1)/N,DIM=ndim-1,NCOPIES=N)
Info%bc(2)%p(1:ny,1:nz)=u(nx+1,:,:)-ut(:,:,2)-&
     SPREAD(SUM(u(nx+1,:,:)-ut(:,:,2),ndim-1)/N,DIM=ndim-1,NCOPIES=N)

IF(ndim==2) THEN
! bottom and top
   Info%bc(3)%p(1:nx,1:nz)=v(1:nx,1,1:nz)-vt(:,:,1)
   Info%bc(4)%p(1:nx,1:nz)=v(1:nx,ny+1,1:nz)-vt(:,:,2)
ELSE
   ! front and back
   Info%bc(3)%p(1:nx,1:nz)=v(:,1,:)-vt(:,:,1)-&
        SPREAD(SUM(v(:,   1,:)-vt(:,:,1),ndim-1)/N,DIM=ndim-1,NCOPIES=N)
   Info%bc(4)%p(1:nx,1:nz)=v(:,ny+1,:)-vt(:,:,2)-&
        SPREAD(SUM(v(:,ny+1,:)-vt(:,:,2),ndim-1)/N,DIM=ndim-1,NCOPIES=N)

! top and bottom
   Info%bc(5)%p(1:nx,1:ny)=w(1:nx,1:ny,1)-wt(:,:,1)
   Info%bc(6)%p(1:nx,1:ny)=w(1:nx,1:ny,nz+1)-wt(:,:,2)

ENDIF 

   CALL take_div(Info,.true.)
!!$IF(ABS(SUM(Info%div))>1.d-8) THEN
!!$   PRINT*,'the averaged divergence is non zero'
!!$   PRINT*,'solution non converging'
!!$   PRINT*,'Sum of divergence is',SUM(Info%div)

!!$ENDIF

div=>Info%div

! modify upper boundary condition
IF(ndim==2) THEN
Info%bc(4)%p(1:nx,1:nz)=Info%bc(4)%p(1:nx,1:nz)+ &
     SUM(div,ndim)*Info%dX(ndim)
!!$     SUM(div,ndim)/N
v(:,ny+1,1)=v(:,ny+1,1)-&
SUM(div(:,:,1),ndim)*Info%dX(ndim)
!!$SUM(div(:,:,1),ndim)/N
ELSE
Info%bc(6)%p(1:nx,1:ny)=Info%bc(6)%p(1:nx,1:ny)+&
     SUM(div,ndim)*Info%dX(ndim)
!!$     SUM(div,ndim)/N
w(:,:,nz+1)=w(:,:,nz+1)-&
     SUM(div,ndim)*Info%dX(ndim)
!!$SUM(div,ndim)/N
END IF
div=div-SPREAD(SUM(div,ndim)/N,DIM=ndim,NCOPIES=N)

!!$
!!$PRINT*,'divergence stat'

!!$PRINT*,'max divergence before = ',maxval(abs(div))
!!$PRINT*,'rms divergence before = ',sqrt(SUM(div**2)/(nx*ny*nz))
!!$PRINT*,'maxval vert sum div = ',maxval(SUM(div,2))

! solve the vertical problems
p=zero

!!$error=POISSON_SOLVERS_THIN(Info,SOLVE,OP=CSR_PRECOND)
!!$STOP ! TO BE REMOVED WHEN EVERYTHING WORKS


error=POISSON_TRIDIAG(Info)





! now we apply the bc to the potential. Note that to generate the off-diagonal
! terms we also need info from pressure at previous level, stored in pt


! apply the operator that adds the bc to the ghost layers, including p from
! previous level
dummy=>adv(Info%G(EXTR_BC)%crsptr,p(1:nx,1:ny,1:nz)+pt(1:nx,1:ny,1:nz),(/2,2,2/))
p=p+dummy; DEALLOCATE(dummy);NULLIFY(DUMMY)
! now subtract from the ghost layers pt
p(0,:,:)=p(0,:,:)-pt(1,:,:)
p(nx+1,:,:)=p(nx+1,:,:)-pt(nx,:,:)
p(:,0,:)=p(:,0,:)-pt(:,1,:)
p(:,ny+1,:)=p(:,ny+1,:)-pt(:,ny,:)
IF(ndim==3) THEN
   p(:,:,0)=p(:,:,0)-pt(:,:,1)
   p(:,:,nz+1)=p(:,:,nz+1)-pt(:,:,nz)
ENDIF
! finally, need to subtract bottom plane to remove noise
IF(ndim==2) THEN
   p(:,:,1)=p(:,:,1)-SPREAD(p(:,1,1),DIM=ndim,NCOPIES=N+2)
ELSE
   p=p-SPREAD(p(:,:,1),DIM=ndim,NCOPIES=N+2)
ENDIF
! now project using this p (.false. since we have already applied bc)
! *****************************
! filter in x
! *****************************


CALL project(Info,.false.)
! save p
pt=p;



! **************************** end vertical *************************

! **************************** begin horizontal *********************
p=zero
! apply averaged boundary conditions
! left and right
Info%bc(1)%p(1:ny,1:nz)=&
     SPREAD(SUM(u(1,:,:)-ut(:,:,1),ndim-1)/N,DIM=ndim-1,NCOPIES=N)

Info%bc(2)%p(1:ny,1:nz)=&
     SPREAD(SUM(u(nx+1,:,:)-ut(:,:,2),ndim-1)/N,DIM=ndim-1,NCOPIES=N)

IF(ndim==2) THEN
! bottom and top
   Info%bc(3)%p(1:nx,1:nz)=v(1:nx,1,1:nz)-vt(:,:,1)
   Info%bc(4)%p(1:nx,1:nz)=v(1:nx,ny+1,1:nz)-vt(:,:,2)
ELSE
   ! front and back
   Info%bc(3)%p(1:nx,1:nz)=&
        SPREAD(SUM(v(:,   1,:)-vt(:,:,1),ndim-1)/N,DIM=ndim-1,NCOPIES=N)
   Info%bc(4)%p(1:nx,1:nz)=&
        SPREAD(SUM(v(:,ny+1,1:nz)-vt(:,:,2),ndim-1)/N,DIM=ndim-1,NCOPIES=N)

! top and bottom
   Info%bc(5)%p(1:nx,1:ny)=w(1:nx,1:ny,1)-wt(:,:,1)
   Info%bc(6)%p(1:nx,1:ny)=w(1:nx,1:ny,nz+1)-wt(:,:,2)

ENDIF
! take divergence
CALL take_div(Info,.true.)
div=>Info%div
!!$PRINT*,'divergence stat'

!!$PRINT*,'max divergence before = ',maxval(abs(Info%div))
!!$PRINT*,'rms divergence before = ',sqrt(SUM(Info%div**2)/(nx*ny*nz))

! average vertically divergence
IF(ndim==2) THEN
   div(:,1,1)=SUM(div(:,:,1),ndim)/N

ELSE
   div(:,:,1)=SUM(div(:,:,:),ndim)/N   
ENDIF
!!$div(2:nx-1,:,:)=(div(1:nx-2,:,:)+2*div(2:nx-1,:,:)+div(3:nx,:,:))/4.d0

! solve D-1 problem
Error=POISSON_SOLVERS_THIN(Info,SOLVE)

! project the field using this pressure
! first add boundary conditions horizontally
IF(ndim==2) THEN
   dummy=>adv(Info%G(10)%crsptr,p(1:nx,1:1,1:1),(/2,2,2/))
   p(:,:,1)=SPREAD(p(:,1,1)+dummy(:,2,2),DIM=2,NCOPIES=ny+2) ! the bounds of dummy are offset relative to p
   DEALLOCATE(dummy);NULLIFY(dummy)
ELSE
   dummy=>adv(Info%G(10)%crsptr,p(1:nx,1:ny,1:1),(/2,2,2/))
   p(:,:,:)=SPREAD(p(:,:,1)+dummy(:,:,2),DIM=3,NCOPIES=nz+2)
   DEALLOCATE(dummy);NULLIFY(dummy)
ENDIF
! **********************
! filter
! **********************
!!$p(1:nx,:,:)=(p(0:nx-1,:,:)+2*p(1:nx,:,:)+p(2:nx+1,:,:))/4.d0

! project 
CALL project(Info,.false.)
! accumulate
pt=pt+p;
! ********************* end horizontal ******************************
pacc=pacc+pt

! check if stop iterate
! this is an expensive step
!!$CALL take_div(Info,.false.)

ITERATE=.true.!sqrt(SUM(Info%div**2)/(nx*ny*nz))>1e-8

i=i+1

!!$PRINT*,' at end of iteration ',i,' l2 norm of div is',sqrt(SUM(div**2)/(nx*ny*nz))

END DO

IF(vs) THEN ! the vertical velocity is calculated integrating the horizontal divergence
      dummy=>adv(Info%G(4)%crsptr,u,(/-1,0,0/))
      IF(ndim==2) THEN

         v(:,1,1)=vt(:,1,1)
         DO j=2,N+1
            v(:,j,1)=v(:,j-1,1)-dummy(:,j-1,1)*Info%dX(ndim)
         ENDDO
         DEALLOCATE(DUMMY);NULLIFY(DUMMY)
      ELSE
         dummy2=>adv(Info%G(5)%crsptr,v,(/0,-1,0/))
         w(:,:,1)=wt(:,:,1)
         DO j=2,N+1
            w(:,:,j)=w(:,:,j-1)-(dummy(:,:,j-1)+dummy2(:,:,j-1))*Info%dX(ndim)
         ENDDO
         DEALLOCATE(dummy,dummy2);Nullify(dummy,dummy2)
      ENDIF
   END IF
   
   Info%p=pacc

!!$CALL take_div(Info,.false.)
!!$div=>Info%div


!!$PRINT*,'final value of Divergence is ',SQRT(SUM(div**2)/(nx*ny*nz))
 END SUBROUTINE hodge_thin

SUBROUTINE hodge_thick(Info,cutoff,ut,vt,wt,solver)
!
TYPE (ELLI):: Info
REAL (KIND=Precision), INTENT (IN) :: ut(1:Info%mX(2),1:Info%mX(3),2)
REAL (KIND=Precision), INTENT (IN) :: vt(1:Info%mX(1),1:Info%mX(3),2)
REAL (KIND=Precision), OPTIONAL, INTENT (IN) :: wt(1:Info%mX(1),1:Info%mX(2),2)
INTEGER, OPTIONAL, INTENT (IN) :: solver
REAL (KIND=Precision), INTENT (IN) :: cutoff
!
LOGICAL :: iterate
INTEGER :: nx,ny,nz,ndim,N,i,err,j
REAL (KIND=Precision), DIMENSION (:,:,:), POINTER :: u,v,w,div,p,dummy
REAL (KIND=Precision), DIMENSION (:,:,:), ALLOCATABLE :: pt
! sanity check
IF(.NOT.PRESENT(wt).AND.Info%ndim==3) THEN
   PRINT*,'wt must be input in 3D'
   STOP
END IF
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3);ndim=Info%ndim;N=Info%mX(ndim)

ALLOCATE(pt(0:nx+1,0:ny+1,0:nz+1))


u=>Info%qe(1:nx+1,1:ny,1:nz,1)
v=>Info%qe(1:nx,1:ny+1,1:nz,2)
w=>Info%qe(1:nx,1:ny,1:nz+1,3)
div=>Info%div
p=>Info%p

Iterate=.true.;i=0
pt=zero
DO WHILE(ITERATE)
   i=i+1


!!$   PRINT*,'++++++++',(SUM(ut(:,1,1))-SUM(ut(:,1,2)))/Info%dx(1)+ &
!!$        (SUM(vt(:,1,1))-SUM(vt(:,1,2)))/Info%dx(2),'++++++++++++'
! set boundary conditions
   Info%bc(1)%p(1:ny,1:nz)=u(1  ,1:ny,1:nz)-ut(:,:,1)
   Info%bc(2)%p(1:ny,1:nz)=u(nx+1,1:ny,1:nz)-ut(:,:,2)

   Info%bc(3)%p(1:nx,1:nz)=v(1:nx,1   ,1:nz)-vt(:,:,1)
   Info%bc(4)%p(1:nx,1:nz)=v(1:nx,ny+1,1:nz)-vt(:,:,2)

   IF(ndim>2) THEN
      Info%bc(5)%p(1:nx,1:ny)=w(1:nx,1:ny,   1)-wt(:,:,1)
      Info%bc(6)%p(1:nx,1:ny)=w(1:nx,1:ny,nz+1)-wt(:,:,2)
   ENDIF

! calulates div (including boundary corrections

   CALL take_div(Info,.true.)
!!$PRINT*, 'max divergence before = ',maxval(abs(Info%div))
!!$
!!$PRINT*,'rms divergence before = ',sqrt(SUM(Info%div**2)/(nx*ny*nz))
!!$PRINT*,'sum of divergence is ',SUM(Info%div)

IF(ABS(SUM(Info%div)*PRODUCT(Info%dX(1:ndim)))>1.d-4) THEN
   PRINT*,'the averaged divergence is non zero'
   PRINT*,'solution non converging'
   PRINT*,'Sum of divergence is',SUM(Info%div),PRODUCT(Info%dX(1:ndim))
   STOP
ENDIF
! solve equation
   IF(PRESENT(solver)) THEN
      err=POISSON_SOLVERS(Info,SOLVE,solver,OP=CSR_PRECOND)
   ELSE
      err=POISSON_SOLVERS(Info,SOLVE,PARDISO_SOLVER)
   ENDIF
   

! subtract gradient
   CALL project(Info,.true.)

   CALL take_div(Info,.false.)
!!$PRINT*, 'max divergence after = ',maxval(abs(Info%div))
!!$
!!$PRINT*,'rms divergence after = ',sqrt(SUM(Info%div**2)/(nx*ny*nz))
!!$PRINT*,'sum of divergence is ',SUM(Info%div)


   ITERATE=sqrt(SUM(Info%div**2)/(nx*ny*nz))>cutoff
   IF(i>10) EXIT
pt=pt+Info%p
END DO
Info%p=pt;
DEALLOCATE(pt)
END SUBROUTINE hodge_thick

FUNCTION POISSON_TRIDIAG(Info)
! Solves the 1D poisson problem d/dz(|g|g^33d/dz)p=div
! it is vectorized in the x (and y) components 
TYPE (Elli) :: Info
INTEGER POISSON_TRIDIAG
!
INTEGER i,j,k,error

SELECT CASE (Info%ndim)

CASE (2)

FORALL (i=1:Info%mX(1))
   Info%p(i,1:Info%mX(2),1)=&
        tridiag(Info%div(i,1:Info%mX(2),1),&
        Info%diag(i,1:Info%mX(2)+1,1),&
        Info%dX(2))
END FORALL

CASE(3)
FORALL (i=1:Info%mX(1),j=1:Info%mX(2))
   Info%p(i,j,1:Info%mX(3))=&
        tridiag(Info%div(i,j,1:Info%mX(3)),&
        Info%diag(i,j,1:Info%mX(3)+1),&
        Info%dX(3))
END FORALL
END SELECT
POISSON_TRIDIAG=ErrOk
END FUNCTION POISSON_TRIDIAG

PURE FUNCTION tridiag(r,f,dz) RESULT (u)
! solves the tridiagonal problem 
! the diagonal of the problem is -a-c where a and c
! are the super and subdiagonal.
! for the problem at hand, a(1:n)=f(1:n)/dz^2
! c(1:n)=f(2:n+1) since f is defined on edges
! also a(1)=0 and c(n)=0 (Neumann conditions)
REAL (KIND=Precision), DIMENSION (:), INTENT (IN) :: r,f
REAL (KIND=Precision), DIMENSION (SIZE(r)) :: u
REAL (KIND=Precision), INTENT (IN) :: dz
!!$INTEGER tridiag
!
INTEGER :: k,n
REAL (KIND=Precision), DIMENSION (size(r)) :: gam
REAL (KIND=Precision) :: z,dz2
!!$tridiag=ErrOk
dz2=dz**2
n=SIZE(u)


   z=-dz2/(f(2))
   u(1)=r(1)*z
   gam(1)=f(2)/dz2*z
   DO k=2,n-1
      z=dz2/(-f(k)*(1+gam(k-1))-f(k+1))
      u(k)=(r(k)-f(k)/dz2*u(k-1))*z
      gam(k)=f(k+1)/dz2*z
   END DO
   ! IF BC are non NEU-NEU, Above loop can be extended to n
   gam(n)=-(2*f(n))/dz2
   z=1./(gam(n)-f(n)/dz2*gam(n-1))
   u(n)=(r(n)-f(n)/dz2*u(n-1))*z
   DO k=n-1,1,-1
      u(k)=u(k)-gam(k)*u(k+1)
   END DO


END FUNCTION tridiag

END MODULE Elliptic_Solvers
